Datos

# Bibliotecas
library(tidyverse)
library(readxl)

# Carga de datos
datos <- read_excel("~/Desktop/Artículos Scientific reports/Resubmission 2/Data/Datos/Clean data/ACP Analysis.xlsx") %>% 
  mutate(Moments = factor(Moments,
                          levels = c("Habituation", "Treatment", 
                                     "Break", "Final")))
datos

Componentes Principales

Retención de variabilidad

library(FactoMineR)
library(factoextra)
library(DT)
acp <- PCA(
  X = datos %>% select(Active:Distressed) %>%
    select(-c(Apathetic, Distressed, Frustrated, Bored)),
  scale.unit = TRUE,
  graph = FALSE
)

as.data.frame(acp$eig) %>% 
  mutate_if(is.numeric, round, digits = 2) %>% 
  rownames_to_column(var = "component") %>% 
  datatable(extensions = 'Buttons', options = list(
     dom = 'Bfrtip',
     buttons = c('csv', 'excel', 'pdf')),
     rownames = FALSE)
  • Gráfico:
fviz_eig(acp, addlabels = TRUE)

Contribución de variables

  • Tabla:
as.data.frame(acp$var$contrib) %>% 
  rownames_to_column(var = "variable") %>% 
  mutate_if(is.numeric, round, digits = 2) %>% 
  datatable(extensions = 'Buttons', options = list(
     dom = 'Bfrtip',
     buttons = c('csv', 'excel', 'pdf')),
     rownames = FALSE)
  • Gráfico:
library(ggsci)
library(tidytext)
as.data.frame(acp$var$contrib) %>% 
  rownames_to_column(var = "adjetivo") %>% 
  pivot_longer(cols = -adjetivo, names_to = "componente",
               values_to = "valor") %>% 
  mutate(componente = as_factor(componente),
         adjetivo = reorder_within(adjetivo, valor, componente)) %>% 
  ggplot(aes(x = adjetivo, y = valor, fill = componente,
             color = componente)) +
  facet_wrap(~componente, scales = "free") +
  geom_col(position = "dodge", alpha = 0.7) + 
  coord_flip() +
  scale_x_reordered() + 
  scale_color_npg() +
  scale_fill_npg() +
  labs(color = "", fill = "") +
  theme(legend.position = "top")

Correlaciones

  • Tabla:
as.data.frame(acp$var$cor) %>% 
  rownames_to_column(var = "variable") %>% 
  mutate_if(is.numeric, round, digits = 2) %>% 
  datatable(extensions = 'Buttons', options = list(
     dom = 'Bfrtip',
     buttons = c('csv', 'excel', 'pdf')),
     rownames = FALSE)
  • Gráfico:
as.data.frame(acp$var$cor) %>% 
  rownames_to_column(var = "variable") %>% 
  mutate(across(where(is.numeric), round, digits = 2)) %>% 
  pivot_longer(cols = -variable) %>% 
  ggplot(aes(x = variable, y = name, fill = value)) +
  geom_tile(color = "black") +
  geom_text(aes(label = value), color = "black") +
  scale_fill_gsea() +
  labs(x = "", y = "", fill = "Correlación:") +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 45, hjust = 1))

Gráficos de componentes

CP1 vs PC2

fviz_pca_biplot(
  acp,
  col.ind = datos$Moments,
  palette = "npg",
  addEllipses = TRUE,
  label = "var",
  col.var = "black",
  repel = TRUE,
  legend.title = "Moments:",
  ellipse.type = "norm",
  ellipse.level = 0.5,
  select.var = list(contrib = 10),
  geom.ind = "point",
  pointshape = 21,
  pointsize = 2.5,
  fill.ind = datos$Moments
)  +
  theme_minimal() +
  theme(legend.position = "top") +
  labs(x = "Component 1", y = "Component 2",
       title = "")

CP1 vs PC3

fviz_pca_biplot(
  acp,
  axes = c(1, 3),
  col.ind = datos$Moments,
  palette = "npg",
  addEllipses = TRUE,
  label = "var",
  col.var = "black",
  repel = TRUE,
  legend.title = "Moments:",
  ellipse.type = "norm",
  ellipse.level = 0.5,
  select.var = list(contrib = 10),
  geom.ind = "point",
  pointshape = 21,
  pointsize = 2.5,
  fill.ind = datos$Moments
)  +
  theme_minimal() +
  theme(legend.position = "top") +
  labs(x = "Component 1", y = "Component 3",
       title = "")

CP2 vs PC3

fviz_pca_biplot(
  acp,
  axes = c(2, 3),
  col.ind = datos$Moments,
  palette = "npg",
  addEllipses = TRUE,
  label = "var",
  col.var = "black",
  repel = TRUE,
  legend.title = "Moments:",
  ellipse.type = "norm",
  ellipse.level = 0.5,
  select.var = list(contrib = 10),
  geom.ind = "point",
  pointshape = 21,
  pointsize = 2.5,
  fill.ind = datos$Moments
)  +
  theme_minimal() +
  theme(legend.position = "top") +
  labs(x = "Component 2", y = "Component 3",
       title = "")

CP1 vs CP2 vs CP3 - Estático

cp1 <- acp$ind$coord[, 1]
cp2 <- acp$ind$coord[, 2]
cp3 <- acp$ind$coord[, 3]

colors <- pal_npg("nrc")(4)
colors <- colors[as.numeric(datos$Moments)]

library(scatterplot3d)
s3d <- scatterplot3d(
  x = cp1,
  y = cp2,
  z = cp3,
  pch = 19,
  color = colors,
  grid = TRUE,
  xlab = "CP1",
  ylab = "CP2",
  zlab = "CP3",
  box = FALSE
)
legend(s3d$xyz.convert(8, 4, 5.5),
       legend = levels(datos$Moments),
       col = pal_npg("nrc")(4),
       pch = 19, box.col = "white")

CP1 vs CP2 vs CP3 - Interactivo

datos$cp1 <- cp1
datos$cp2 <- cp2
datos$cp3 <- cp3

library(plotly)
fig <- plot_ly(data = datos,
               x = ~cp1, y = ~cp2, z = ~cp3,
               color = ~Moments,
               colors = pal_npg("nrc")(4))
fig <- fig %>% add_markers()
fig

Inferencia

Distribuciones

datos %>% 
  select(Moments, cp1:cp3) %>% 
  pivot_longer(cols = -Moments) %>% 
  ggplot(aes(x = value)) +
  facet_wrap(~name) +
  geom_density()

Boxplots

datos %>% 
  select(Moments, cp1:cp3) %>% 
  pivot_longer(cols = -Moments) %>% 
  ggplot(aes(x = Moments, y = value, color = name, fill = name)) +
  geom_boxplot(alpha = 0.5) +
  scale_fill_npg() +
  scale_color_npg() +
  labs(color = "", fill = "")

Anova CP1

  • Resultado análisis de varianza:
anova_cp1 <- aov(cp1 ~ Moments, data = datos)
tidy(anova_cp1)
  • Estimación de medias:
library(emmeans)
emmeans(anova_cp1, specs = "Moments")
##  Moments     emmean    SE  df lower.CL upper.CL
##  Habituation -0.989 0.486 281   -1.945  -0.0333
##  Treatment    1.164 0.209 281    0.753   1.5757
##  Break       -0.846 0.243 281   -1.324  -0.3684
##  Final       -1.913 0.486 281   -2.869  -0.9574
## 
## Confidence level used: 0.95
  • Diferencias entre medias - Tabla:
tidy(TukeyHSD(anova_cp1))
  • Diferencias entre medias - Gráfico:
tidy(TukeyHSD(anova_cp1)) %>% 
  ggplot(aes(x = contrast, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_point() +
  geom_errorbar(width = 0.2) +
  geom_hline(yintercept = 0, linetype = 2, color = "red") +
  coord_flip()

  • Residuales:
par(mfrow = c(2, 2))
plot(anova_cp1)

Anova CP2

  • Resultado análisis de varianza:
anova_cp2 <- aov(cp2 ~ Moments, data = datos)
tidy(anova_cp2)
  • Estimación de medias:
emmeans(anova_cp2, specs = "Moments")
##  Moments     emmean    SE  df lower.CL upper.CL
##  Habituation -1.009 0.331 281   -1.660   -0.358
##  Treatment    0.908 0.142 281    0.628    1.188
##  Break       -0.657 0.165 281   -0.983   -0.332
##  Final       -1.264 0.331 281   -1.915   -0.613
## 
## Confidence level used: 0.95
  • Diferencias entre medias - Tabla:
tidy(TukeyHSD(anova_cp2))
  • Diferencias entre medias - Gráfico:
tidy(TukeyHSD(anova_cp2)) %>% 
  ggplot(aes(x = contrast, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_point() +
  geom_errorbar(width = 0.2) +
  geom_hline(yintercept = 0, linetype = 2, color = "red") +
  coord_flip()

  • Residuales:
par(mfrow = c(2, 2))
plot(anova_cp2)

Anova CP3

  • Resultado análisis de varianza:
anova_cp3 <- aov(cp3 ~ Moments, data = datos)
tidy(anova_cp3)
  • Estimación de medias:
emmeans(anova_cp3, specs = "Moments")
##  Moments      emmean    SE  df lower.CL upper.CL
##  Habituation -0.2515 0.240 281  -0.7231   0.2201
##  Treatment    0.1786 0.103 281  -0.0243   0.3816
##  Break       -0.1892 0.120 281  -0.4250   0.0466
##  Final        0.0438 0.240 281  -0.4278   0.5154
## 
## Confidence level used: 0.95
  • Diferencias entre medias - Tabla:
tidy(TukeyHSD(anova_cp3))
  • Diferencias entre medias - Gráfico:
tidy(TukeyHSD(anova_cp3)) %>% 
  ggplot(aes(x = contrast, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_point() +
  geom_errorbar(width = 0.2) +
  geom_hline(yintercept = 0, linetype = 2, color = "red") +
  coord_flip()

  • Residuales:
par(mfrow = c(2, 2))
plot(anova_cp3)

Modelos Mixtos

Variación individual (1)

  • CP1:
datos %>%  
  mutate(Individuo = as_factor(Individuo)) %>% 
  group_by(Moments, Individuo) %>% 
  summarise(promedio_cp = mean(cp1)) %>% 
  ggplot(aes(x = Moments, y = promedio_cp, color = Individuo, group = Individuo)) +
  geom_point() +
  geom_line() +
  labs(color = "") +
  guides(colour = guide_legend(nrow = 3)) +
  theme(legend.position = "top")

  • CP2:
datos %>%  
  mutate(Individuo = as_factor(Individuo)) %>% 
  group_by(Moments, Individuo) %>% 
  summarise(promedio_cp = mean(cp2)) %>% 
  ggplot(aes(x = Moments, y = promedio_cp, color = Individuo, group = Individuo)) +
  geom_point() +
  geom_line() +
  labs(color = "") +
  guides(colour = guide_legend(nrow = 3)) +
  theme(legend.position = "top")

Variación individual (2)

  • CP1:
datos %>%  
  mutate(Individuo = as_factor(Individuo)) %>% 
  group_by(Moments, Individuo) %>% 
  summarise(promedio_cp = mean(cp1)) %>% 
  ggplot(aes(x = Moments, y = promedio_cp, color = Individuo, group = Individuo)) +
  facet_wrap(~Individuo, scales = "free", ncol = 5) +
  geom_point() +
  geom_line() +
  labs(color = "") +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1))

  • CP2:
datos %>%  
  mutate(Individuo = as_factor(Individuo)) %>% 
  group_by(Moments, Individuo) %>% 
  summarise(promedio_cp = mean(cp2)) %>% 
  ggplot(aes(x = Moments, y = promedio_cp, color = Individuo, group = Individuo)) +
  facet_wrap(~Individuo, scales = "free", ncol = 5) +
  geom_point() +
  geom_line() +
  labs(color = "") +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1))

datos %>% filter(Pieza == "Pieza4")

Modelo CP1

  • Análisis de varianza:
library(nlme)
mod_cp1 <- lme(cp1 ~ Moments, method = "REML", 
               random = ~ Moments | Individuo,
               data = datos,
               control = lmeControl(opt = "optim"))

tidy(car::Anova(mod_cp1, type = "II"))
  • Estimación de medias:
emmeans(mod_cp1, specs = "Moments")
##  Moments     emmean    SE df lower.CL upper.CL
##  Habituation  -1.06 0.453 29   -1.989   -0.135
##  Treatment     1.22 0.313 29    0.582    1.864
##  Break        -0.97 0.232 29   -1.446   -0.495
##  Final        -1.99 0.453 29   -2.913   -1.059
## 
## Degrees-of-freedom method: containment 
## Confidence level used: 0.95
  • Diferencias entre medias - Tabla:
tidy(pairs(emmeans(mod_cp1, ~Moments))) %>% 
  dplyr::select(contrast, estimate, std.error, adj.p.value)
  • Diferencias entre medias - Gráfico:
library(multcomp)
library(multcompView)

medias <- as.data.frame(cld(emmeans(mod_cp1, ~Moments), 
                            alpha = 0.05, Letters = letters,
                            adjust = "tukey"))


medias  %>% 
  ggplot(aes(x = Moments, y = emmean, ymin = emmean - SE,
             ymax = emmean + SE)) +
  geom_point(position = position_dodge(0.9), alpha = 1) +
  geom_errorbar(width = 0.1) +
  geom_line(aes(group = 1)) +
  geom_text(aes(label = .group, y = emmean),
            position = position_dodge(0.9),
            show.legend = FALSE, color = "black",
            hjust = -0.3,
            vjust = -0.5) +
  scale_color_grey() +
  scale_fill_grey() +
  geom_hline(yintercept = 0, color = "black", linetype = 2) +
  labs(y = "", x = "Moments") +
  theme_bw() +
  theme(legend.position = "none") 

  • Residuales:
par(mfrow = c(1, 2))
#plot(mod_cp1)

residuales <- mod_cp1$residuals[, 1]
ajustados <- mod_cp1$fitted[, 1]

qqnorm(residuales)
qqline(residuales)

plot(x = ajustados, y = residuales)
abline(h = 0, col = "red", lty = 2)

Modelo CP2

  • Análisis de varianza:
mod_cp2 <- lme(cp2 ~ Moments, method = "REML", 
               random = ~ Moments | Individuo,
               data = datos,
               control = lmeControl(opt = "optim"))

tidy(car::Anova(mod_cp2, type = "II"))
  • Estimación de medias:
emmeans(mod_cp2, specs = "Moments")
##  Moments     emmean    SE df lower.CL upper.CL
##  Habituation -1.009 0.316 29   -1.657   -0.362
##  Treatment    0.967 0.185 29    0.588    1.347
##  Break       -0.657 0.158 29   -0.981   -0.334
##  Final       -1.288 0.323 29   -1.948   -0.627
## 
## Degrees-of-freedom method: containment 
## Confidence level used: 0.95
  • Diferencias entre medias - Tabla:
tidy(pairs(emmeans(mod_cp2, ~Moments))) %>% 
  dplyr::select(contrast, estimate, std.error, adj.p.value)
  • Diferencias entre medias - Gráfico:
medias2 <- as.data.frame(cld(emmeans(mod_cp2, ~Moments), 
                            alpha = 0.05, Letters = letters,
                            adjust = "tukey"))


medias2  %>%  
  ggplot(aes(x = Moments, y = emmean, ymin = emmean - SE,
             ymax = emmean + SE)) +
  geom_point(position = position_dodge(0.9), alpha = 1) +
  geom_errorbar(width = 0.1) +
  geom_line(aes(group = 1)) +
  geom_text(aes(label = .group, y = emmean),
            position = position_dodge(0.9),
            show.legend = FALSE, color = "black",
            hjust = -0.3,
            vjust = -0.5) +
  scale_color_grey() +
  scale_fill_grey() +
  geom_hline(yintercept = 0, color = "black", linetype = 2) +
  labs(y = "", x = "Moments") +
  theme_bw() +
  theme(legend.position = "none") 

  • Residuales:
par(mfrow = c(1, 2))
#plot(mod_cp1)

residuales2 <- mod_cp2$residuals[, 1]
ajustados2 <- mod_cp2$fitted[, 1]

qqnorm(residuales2)
qqline(residuales2)

plot(x = ajustados2, y = residuales2)
abline(h = 0, col = "red", lty = 2)

Exportando gráficos

# Exportando gráfico pdf
ggsave(filename = "graphics_paper/acp.pdf", 
       plot = g1,
       device = "pdf",
       dpi = 96,
       units = "cm",
       width = 35.14, # igual a 531px
       height = 14.05) # igual a 1328px

# Exportando gráfico Tiff
ggsave(filename = "graphics_paper/acp.tiff", 
       plot = g1,
       device = "tiff",
       dpi = 96,
       units = "cm",
       width = 35.14, # igual a 531px
       height = 14.05) # igual a 1328px